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We carry out a perturbative analysis, of high order in the tunneling parameter, of the ground state 
of the homogeneous Bose-Hubbard model in the Mott insulator phase. This is made possible by a 
diagrammatic process chain approach, derived from Kato's representation of the many-body per- 
turbation series, which can be implemented numerically in a straightforward manner. We compute 
ground-state energies, atom-atom correlation functions, density-density correlations, and occupation 
number fluctuations, for one-, two-, and three-dimensional lattices with arbitrary integer filling. A 
phenomenological scaling behavior is found which renders the data almost independent of the filling 
factor. In addition, the process chain approach is employed for calculating the boundary between 
the Mott insulator phase and the superfluid phase with high accuracy. We also consider systems 
with dimensionalities d > 3, thus monitoring the approach to the mean-field limit. The versatility 
of the method suggests further applications to other systems which are less well understood. 
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I. INTRODUCTION 

With the seminal experiment by Greiner et al.j- who 
have observed the quantum phase transition from a Mott 
insulator to a superfluid^ in a gas of ultracold Rubidium 
atoms stored in a three-dimensional optical lattice, the 
regime of strongly interacting Bose gases has become ac- 
cessible. Due to its perfectly controllable parameters, 
the experimental set-up, as suggested by Jaksch et al.,^ 
provides a fantastic testing ground for quantum many- 
body physicsi^ Meanwhile, the transition has also been 
observed in one- and two-dimensional optical lattices.^'- 
Ultracold atoms in optical lattices are described by the 
Bose-Hubbard model, ^ which incorporates two compet- 
ing trends in an elementary manner: While the repulsive 
interaction between the atoms tends to localize the par- 
ticles at individual sites of the lattice potential, tunnel- 
ing between neighboring sites favors delocalization, and 
tends to suppress phase fluctuations. 

The Bose-Hubbard model and its descendants have 
been intensively studied within the last years. Important 
techniques for monitoring its ground-state properties and 
the phase diagram include the mean-fleld approach;^ 
strong-coupling expansions,.^ methods using the density 
matrix renormalization group (DMRG) )^°i"'^^'^^i^^ and, 
more recently, quantum Monte Carlo (QMC) simula- 
tions 

In the present work we apply a recently suggested pro- 
cess chain approachlS^ to the d-dimensional homogeneous 
Bose-Hubbard model, in order to investigate the proper- 
ties of its ground state in the Mott insulator phase for 
any integer filling factor g > 1, as well as the phase di- 
agram for d > 2, by means of a high-order expansion 
in the tunneling parameter. To achieve this, we have 



turned Kato's representation of the perturbation series^^ 
into a numerically executable algorithm which handles 
symbolic diagrams as inputs. Order by order, each ob- 
servable then is represented by a set of such diagrams, 
equipped with appropriate weight factors. 

The paper is organized as follows: In Sec.|TT]the Bose- 
Hubbard model is briefly recapitulated. Our zeroth-order 
Hamiltonian contains the local on-site interaction only, 
while the tunneling term will be treated perturbatively. 
Section IIIII introduces Kato's formulation of the pertur- 
bation series, and explains the required elements of the 
diagrammatic process chain approachii^ As an instruc- 
tive example for this technique, we outline in some de- 
tail how to calculate the fourth-order energy correction 
induced by the tunneling term. We also discuss how 
fully correlated ground-state expectation values are ac- 
cessible via the perturbational approach. In Sec. IIVI we 
present results for the ground-state energy, correlation 
functions, and occupation number fluctuations, as ob- 
tained for homogeneous hypercubic lattices with dimen- 
sionality d — 1, 2, and 3. The required diagrams are 
developed, and a phenomenological scaling is suggested, 
which makes the data almost independent of the filling 
factor g. Results typically are calculated up to tenth or- 
der in the tunneling parameter. For d = 1 and unit fill- 
ing we obtain perfect agreement with the results of the 
high-order symbolic perturbative expansion of Damski 
and Zakrzewski.-*^ We follow their instructive work to 
some extent in spirit, opening up the regimes of higher d 
for any g. The Mott insulator-to-superfluid phase transi- 
tion is then discussed in Sec. |V] within the framework of 
the process chain approach. The phase boundary is de- 
termined by invoking the method of the effective poten- 
tial which provides a convenient indicator for the tran- 
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sition point) ^^'^^ this leads to a computational scheme 
which again can be expressed in terms of diagrams^ 
That scheme is worked out exactly in the case of infinite 
lattice dimensionality d, and then reproduces the well- 
known mean- field result. For d < oo our data, obtained 
by numerically evaluating the diagrams up to some or- 
der I/, lend themselves to a simple extrapolation proce- 
dure for v oo. Critical parameters are computed in 
this manner for c? = 2 and d = 3 and any filling factor g; 
for g = 1, they compare quite favorably to recent high- 
precision QMC data.^Siii The last part of this Sec. |V] 
details how the mean-field prediction is approached with 
increasing dimensionality of the system. We close the 
paper with some concluding remarks and a short outlook 
in Sec Ell 

II. THE MODEL 

Ultracold atoms in optical lattices are fairly well de- 
scribed by a single-band Bose-Hubbard model. Ultracold 
in this context means that the thermal de Broglie wave- 
length is at least equal to the lattice constant, i.e., to half 
the wavelength of the laser radiation generating the lat- 
tice. The assumptions underlying this model, requiring 
that the thermal and the interaction energies be smaller 
than the gap between the lowest and the first excited 
Bloch band, are fulfilled if the lattice is sufficiently deep. 
Denoting the interaction energy of a pair of particles oc- 
cupying the same lattice site by U , the chemical potential 
at site i by /ii, and the hopping matrix element connect- 
ing well i and well j by , the model takes the form 

i i,j i 

Here a| and are the bosonic creation and annihilation 
operators for site i, and fii = dla^ is the correspond- 
ing number operator. The chemical potential fii can in- 
corporate an arbitrary confining potential, and then de- 
pends on the lattice site. By choosing appropriate hop- 
ping elements Jij, longer-range or anisotropic hopping 
can be modeled. In this study we stick to the paradig- 
matically simple case with site-independent chemical po- 
tential fj, and isotropic nearest-neighbor hopping of pos- 
itive strength J on a hypercubic lattice. Utilizing the 
interaction energy U as the energy scale of reference, the 
dimensionless Hamiltonian then reads as 

HbH = Hf)+ Htun , (2) 

where 

Ho = ^^ni{hi - 1) - fi/U^Ut (3) 

i i 

is site-diagonal, and 

Htun = ~J/U J2 ajaj (4) 



describes tunneling between adjacent sites, with in- 
dicating that the sum is restricted to nearest neighbors. 
One can easily oversee two limiting cases: If the lattice 
is very deep, tunneling even between neighboring wells is 
inhibited, because the tunneling parameter J/U vanishes 
exponentially with increasing lattice depthi2ii2^ The sites 
then decouple and the Hamiltonian becomes local; one 
only needs to consider Hq. Minimizing the on-site en- 
ergy Si = ni{ni — l)/2 — riifi/U of a system with J/U — 
leads to the site-occupancies^ 

r for ^,/U < 
' \g for g-Kfi/U <g , l<gen . 

Thus, as long as the chemical potential fi/U falls between 
g — 1 and g, each site is occupied by an integer number 
g — N/M of atoms, where N is the total number of 
particles, and M the number of lattice sites. Denoting 
the vacuum state by |0), the Hq ground state then is 
given by the product state 

M (diY 

i-)=n^io)- (6) 

i=l ^ " 

For integer /i/J7 the ground state is 2*^-fold degenerate. 
The parameter regime gJ/U ^ 1 gives rise to insulating 
phases, because the system remains incompressible for 
small, but non-zero tunneling strength. A small change 
of the chemical potential fj, then does not lead to a change 
of the occupation number: d{ni)/d^ = 0. 

The opposite limiting case appears when the interac- 
tion between the particles can be neglected in compari- 
son with the kinetic energy, gJ/U ^ 1. The ground state 
then becomes an ideal Bosc-Einstcin condensate with all 
particles occupying the zero-quasimomentum state of the 
lowest band: 

1 / 1 " \" 

Observe that this is an eigenstate of i/tun- Nonetheless, 
in what follows we use the site-diagonal Hamiltonian Hq 
and the Fock state ([6]) as the starting point for our per- 
turbative analysis. 



III. KATO FORMALISM AND PROCESS 
CHAIN APPROACH 

For calculating corrections to the ground-state energy, 
and further ground-state expectation values, we resort 
to the representation of the perturbation series given by 
K&toJ^^ Its distinct advantage lies in the fact that it 
yields closed expressions for the perturbative corrections 
in any order, in contrast to the more familiar recursive 
formulation of the Rayleigh-Schrodinger perturbation se- 
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The ground state |m) of the Hamihonian Hq is a prod- 
uct of local Fock states with g particles sitting at each 
site. When this system is subjected to some perturba- 
tion V , not necessarily given by i?tun, the nth-order cor- 
rection to its energy is given by the trac o^^i^^ 



tr 



{Ui} 



(8) 



where the sum runs over all possible sequences {ag} of 
n+1 non-negative integers which obey the condition 



n+l 



(9) 



The operators 5" linking the individual perturbation 
events V are given by 



\m) (m| 



(0)^a 



for a = 

for a > (10) 



with the energies E^^ = M [g{g - l)/2 - {^/U)g] and 

Ei^^ — J2i ~ l)/2 — i^J'/U)ni] of the unperturbed 

Hq eigenstates |m) and respectively. Because these 
Fock states form an orthonormal basis, one easily derives 



5° 5" =0 for a > 



(11) 



To see how Eq. ([5]) works, let us consider the energy 
correction in second order, i.e., for n = 2. The partition 
problem ^ then has the solutions {1, 0, 0}, {0, 1, 0} and 
{0,0,1}. Accordingly, one finds 



tr [SWS"VS" + S"VSWS" + S"VS"VS^] 

tr [S°VSWS°] 

{m\VS^V\ni) 

ni\V\i){i\V\ni) 



E 



E, 



(0) 



E.. 



(0) 



(12) 



In the second step, cyclic interchangeability of operators 
under a trace has been used, together with Eq. pT|). The 
final expression p2)) coincides exactly with the famil- 
iar textbook result provided by the Rayleigh-Schrodinger 
perturbation theory (see, e.g., Refs. 26,27), as it should. 

Due to the restriction ([9]), in any order n at least two 
superscripts a£ are equal to zero, so that the trace ([S]) 
can always be rewritten as a sum of matrix elements of 
the standard form {m\VS°'^V . . . VS°'"-^V\m). Such el- 
ements will be called Kato terms. We regard each Kato 
term as a sum over process chains^ leading from |m) 
back to |m), with the individual processes corresponding 
to non-zero matrix elements 



a) 



b) 



FIG. 1: Second-order tunneling processes on a lattice. While 
path a) contributes, path b) gives no contribution to the en- 
ergy correction, as final and initial state do not coincide. 



In particular, when calculating energy corrections we 
identify each process with a term of -fftun, setting 



V 



a] a. 



(13) 



Each Kato term now can be viewed as a sum over cer- 
tain chains of tunneling processes on the lattice. Because 
these Kato terms represent expectation values with re- 
spect to the state |m), each chain has to start and to end 
in the state |m). Thus, only closed loops of tunneling 
processes contribute to the energy correction. A simple 
example may illustrate this fact: Consider the energy 
correction in second order, given by Eq. p^ . One then 
has two tunneling processes, which could take place any- 
where on the lattice. But only process chains for which 
initial and final state both coincide with |m) give a con- 
tribution. This requires to tunnel back and forth, thus 
producing the only closed loop with two individual tun- 
neling processes, as illustrated in Fig. [T] Such closed 
loops of tunneling processes will be denoted as paths in 
the following. The respective sequence of the individual 
processes, i.e., their ordering, is quite important for the 
evaluation of the matrix elements, as will become evident 



The particular perturbation (|13p gives no contributions 
in odd orders, because no closed loops can be formed with 
an odd number of tunneling processes on a cubic lattice. 
In fourth order, the general Kato terms are 

E^J^^ = {m\VS'^VS^VS°V\m) 

+ 2 {m\VS'^VS^VS°V\m) 

+ {m\VS'^VS'^VS^V\m) 

+ {m\VSWS^VS^V\m) . (14) 

The first term requires that the initial state |m) be re- 
covered after the first process (numbered from right to 
left), because = — |m)(m| occurs. With the pertur- 
bation (I13p this is impossible. Hence, when treating per- 
turbations that vanish in first order, like the tunneling 
events (jl3p . one can further reduce the number of Kato 
terms. To the third term in Eq. (|14p only chains revisit- 
ing |m) after two processes contribute, while evaluating 
the fourth term requires to take into account all those 
permutations of the processes forming the closed loop 
which do not feature the state |m) as an intermediate 
state. With four tunneling processes one can form lots of 
closed loops on the lattice, but as the system is homoge- 
neous, paths which are topologically identical contribute 
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a) b) c) 

FIG. 2: Examples of topologically identical paths occurring 
in fourth order when calculating the energy correction. These 
paths are described by the same diagram. 

in the same way. We subsume those topologically iden- 
tical paths under a diagram. Examples of topologically 
identical paths are sketched in Fig. [51 According to the 
linked cluster theoremj2^ disconnected diagrams do not 
contribute. With every diagram we associate a weight 
factor which incorporates, on the one hand, the subsum- 
mation of topologically identical paths and, on the other, 
an additional factor of for a diagram visiting s lattice 
sites, as each of these s sites can be the "origin" of the 
diagram. 

Thus, when calculating energy corrections the num- 
ber V of tunneling processes agrees with the respective 
order n of the perturbation series, and the overall pro- 
gram for determining these corrections on a hypercubic 
lattice to a given order i' consists of the following steps: 

• Generate the Kato terms provided by the pertur- 
bation series ([5]) in j/th order. This step is indepen- 
dent of the particular problem under study: Once 
the Kato terms are known, they can be used for 
all kinds of perturbative calculations. Group these 
terms together as far as possible, taking into ac- 
count that odd orders never contribute, as there 
are no closed loops with an odd number of tunnel- 
ing processes. 

• Create all paths representing a closed loop with v 
tunneling processes. Subsume topologically iden- 
tical paths to diagrams, and append the correct 
weight factors. 

• For each diagram, go through all permutations of 
the individual processes; for each particular se- 
quence thus obtained, determine those Kato terms 
which match it. Compute the corresponding matrix 
elements, including the respective energy denomi- 
nators. Sum up the contributions of all sequences 
and all diagrams. 

In high orders this procedure becomes more and more 
cumbersome, as the order v enters factorially, and both 
the number of diagrams and the number of Kato terms 
grows roughly exponentially with Table [T] demon- 
strates this for dimensionalities d — 1,2, and 3. Nonethe- 
less, a considerable advantage offered by this scheme con- 
sists in the fact that its implementation on a computer 
is straightforward, while a representation of the entire 
Hilbert space is not necessary. The strategy of determin- 
ing the contributions to the perturbation series from all 
possible paths on the lattice then allows us to treat filling 



TABLE I: Number of diagrams required by the energy correc- 
tion for lattice dimensionality d, and number of Kato terms, 
vs. the order v. If the perturbation vanishes to first order, one 
is left with the reduced number of Kato terms. This number 
is further diminished in case of the energy correction (last col- 
umn), because an even number of tunneling processes has to 
appear between two projection operators S". 



1/ 


No. of diagrams 


No. of Kato terms 




d= 1 


d = 2 


d = 3 


general 


reduced 


energy relevant 
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a) b) c) 

FIG. 3: Diagrams contributing in fourth order to the energy 
correction. We focus on diagram a) in our example calcu- 
lation. For lattice dimensionality d, the respective weight 
factors are a) 3d(2d - l)/3, b) 2d/2, and c) 2d(2d - 2)/4. 

factors and dimensionalities which are difficult to reach 
by other approaches. In order to clarify the above steps, 
we give an explicit example. 

A. Example 

Let us determine the fourth-order correction of the 
ground-state energy due to the perturbation given by 
Htun- As this perturbation is not diagonal in the 
Fock basis, the first and the second term in Eq. van- 
ish. The remaining Kato terms are 

eI^'> = {m\VS^VS°VS^V\m) 

+ {m\VSWsWS^V\m) . (15) 

In general, if the perturbation does not contribute to 
first order, the number of Kato terms can be signifi- 
cantly reduced, which leads to a substantial computa- 
tional speedup. For calculating the energy correction 
one can actually reduce the number of terms still further, 
since only an even number of tunneling processes can ap- 
pear between two projection operators 5"°. As shown in 
Tab. m this approximately halves the number of terms 
required in tenth order. In our example, we have to 
evaluate the three diagrams depicted in Fig. [3l and to 
determine their weight factors. 

We restrict ourselves to the calculation of the dia- 
gram a) listed in Fig. [31 and denote the individual tunnel- 
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FIG. 4: Fourth-order diagram with individual tunneling pro- 
cesses labeled a, b, c, and d, as considered in the example 
calculation. 



In the same manner, the other 22 permutations of the 
processes a, b, c, and d have to be evaluated. Summing up 
all the resulting contributions, multiplying by the weight 
factor d(2d—l) pertaining to this particular diagram, and 
then adding the other two diagrams with their respective 
weight factors, one arrives at the total energy correction 
in fourth order. 



ing processes by a, 6, c, and d, as indicated in Fig.[4l For 
this computation a system with three lattice sites suffices. 
The ground state then is represented by |m) — \g,g,g), 
with filling factor g. Out of a total number of 4! = 24 
permutations, the first sequence to be treated here is 
{a,b,c,d), leading to the following succession of inter- 
mediate states: 

\9,9,9) ^ |5 - 1,.9+ 1,5) 

^15-1,5+1,5)^15,5,5)- (16) 
Invoking the familiar ladder relations 
d \n) = \fn\n — 1) 

a^\n) = x/irTT|n+ 1) (17) 

for bosonic annihilation and creation operators, the fac- 
tors acquired by tunneling combine to gig -I- 1)'^(J/C/)^. 
Because the initial state does not occur as an intermedi- 
ate state here, this particular sequence does not match 
the first term in Eq. so that only the second one con- 
tributes. As one particle- hole pair is present in each inter- 
mediate state, the three individual energy denominators 
are —Ef^ = — 1 (in multiples of the pair- interaction 
energy [/). The full energy denominator therefore is 
given by (— 1)(— 1)(— 1) = —1. Thus, the contribution to 
the energy correction provided by the sequence (a, 5, c, c?) 
is 



A£;(a^b,c,d) - -5(5 + 1)^ 



(18) 



Next, we treat the sequence (o, c?, 6, c): 
15,5,5) ^ I5 - 1,5+1,5) 



1,9,9) 



^ |5,5- 1,5+ 1) ^ 15,5,5) 



(19) 



Here the initial state is recovered after the second tun- 
neling process, leading to a contribution of the first term 
of Eq. (fT5|) , whereas the second one does not match. The 
prefactor due to tunneling now is g^(g + 1)^(J/J7)^; the 
energy denominator becomes (— 1)^(— 1) = —1. Since S° 
yields another factor of —1 (see Eq. (fTO)) ). the contribu- 
tion of this sequence reads as 



A£^(a,<i,6,c) = 5^(5- 



(20) 



B. Ground-state expectation values 

The technique introduced above also allows one to cal- 
culate expectation values (Hi) of observables Hi in the 
ground state of the homogeneous Bose-Hubbard model as 
expansions in the tunneling strength J/U. Considering 
an extended Hamiltonian 



H — Hq + XHt, 



(21) 



its ground-state energy generically possesses an expan- 
sion of the form 



(n,m) 



(22) 



Denoting the ground-state wave function of the full sys- 
tem (j2ip by |'!/;(A, 77)), the Hellmann-Feynman theorem 
states 

-E={i;{X,rj)\ — mX,rj)), (23) 



implying 



mA"ry™-i£;("^") = (^j{X,r])\Hi\^j{X,'n)) (24) 



and thus resulting in 



= (V(i, 0)1^11^^(1,0)) EE (Fi 



(25) 



This means that an implementation of Kato's perturba- 
tion series can be used for computing the desired ground- 
state expectation values by considering the perturbation 
V — Htun + Hi to first order in Hi : Each process chain 
appearing to order n = v + I m the perturbation series 
then contains v tunneling events described by -fftun, and 
only one process Hi. 

For further details concerning the process chain ap- 
proach we refer to Ref . [H. 



IV. RESULTS 
A. Energy corrections 

As the preceding example has shown, the process chain 
approach in principle works for lattices of any dimension- 
ality, with arbitrary filling factor g, but with increasing 
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FIG. 5: (Color online) Logarithm of the coefRcients a'"' 
for the 3D Bose-Hubbard system with filling factors g = 
1, 2, 3, 5, 10, 20, 50, 100, as defined in Eq. (ge)). The inset shows 
the logarithm of the scaled coefficients (|27p , which are almost 
independent of g. The coefficients grow approximately expo- 
nentially with the order v. The results for the 2D system are 
qualitatively similar to these. 
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FIG. 6: (Color online) Ground-state energy correction per site 
{E-E^^)/M in multiples of the pair interaction energy U for 
filling factors g = 1, 2, 5, 10, 20, 50, 100 for the 3D (upper data 
points), 2D (lower data points), and the ID system (inset). 
Vertical lines indicate the respective critical value {J/U)c for 
the Mott insulator-to-superfluid transition with g = 1 (see 
SeclVjl. As a result of the scaling (I29p . data points for differ- 
ent g fall almost on top of each other. 



order it quickly becomes impracticable to write down the 
resulting terms by hand, since their number proliferates 
rapidly, and it is unlikely that they combine to yield a 
simple expression. However, a numerical implementation 
on a computer is technically feasible and straightforward. 

With our current implementation we are able to calcu- 
late energy corrections per lattice site for the one-, two-, 
and three-dimensional (ID, 2D, and 3D) Bose-Hubbard 
model up to 10th order (12th order in the ID case) in the 
tunneling coupling J/U for any integer filling factor g in 
the form 



E-E, 



(0) 



u=2 



(26) 



Our data for the ID system with g — I agree accurately 
with the results reported by Damski and Zakrzewski)2£ 
who have performed a high-order symbolic perturbative 
expansion for this particular situation. In the 2D and 
the 3D case the coefRcients a*-''^(g) grow to good approx- 
imation exponentially with the order v, as Fig. O demon- 
strates. It is of interest to observe that scaling these 
coefRcients a'^'^\g) by factors g{g + 1) leads to data 



aM(.g) 



V5(5+l) 



(27) 



which are almost independent of the filling factor, as wit- 
nessed by the inset in Fig. [5l This is intuitively intelli- 
gible, since \/3(5 + 1) is a typical factor accompanying 
a tunneling process on a lattice which contains g parti- 
cles per site on the average. Therefore, one can approx- 
imately transform the coefRcients a'^'^\gi) pertaining to 
one Riling factor gi to those referring to another factor 172: 



(51) 



.91(51 



1) 



52(52 + 1) 



u/2 



(52) . (28) 



This relation even is exact in second order, whereas small 
deviations occur in higher orders, for which it still re- 
mains a very good estimate. Naturally, the largest devi- 
ation from this scaling behavior occurs for Riling factor 
g = 1, as will repeatedly become visible in our data. 
With this observation in mind, wc write 



E-E, 
W 



(0) 



= -E«^'^M ^5(5 + 1 



i/=2 



(29) 



Hence, when plotting in Fig. [H] the energy corrections 
as functions of \J g{g + 1)J/U, graphs originating from 
different Riling factors g practically lie on to p of each 
other. Of course, when keeping y^g{g + I) J/U con- 
stant while increasing g, the zeroth-order term Em /M — 
g{g— l)/2 — ^/U becomes dominant, and the corrections 
become relatively small. 



Atom-atom correlation function 



The atom-atom correlation function is deRned by 



C(f* j) = = {ajciA 



(30) 



with a lattice vector f^j- pointing from site j to site i. 
Since the system is invariant under translations by in- 
teger multiples of lattice vectors, Cij depends only on 
Tij , but not on the individual sites i and j. In the limit 
J/U the correlation function Cj j vanishes for any 
non-zero r,;.j , whereas one has Cij = g in the BEG limit 
J/U 00, when all particles condense into the low- 
est Bloch state.?S The calculation of the expectation val- 
ues (1301) provides an example of the strategy outlined in 
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a) b) c) 

FIG. 7: Diagrams of fourth order in J/17 required for calculat- 
ing C{'r) with r = [1, 1,0]. The associated weight factors are 
a) 12d — 20, b) 12d — 8, and c) 4. The dashed arrow, pointing 
from site j to site i, describes the action of the operator (|31|l . 
while the solid arrows again correspond to nearest-neighbor 
tunneling processes Q. 

Subsec. IIIIBl setting 

H, = ala^ . (31) 

Hence, we can employ the same implementation of the 
perturbation series as used for the energy correction; only 
the diagrams have to be adapted. Each process chain now 
has to contain one process (pij) . which has to be strictly 
distinguished from the nearest-neighbor tunneling pro- 
cesses described by -fftun even if i and j label adjacent 
sites. We therefore depict this process Hi by a dashed 
arrow. In Fig. [7] we display the diagrams contributing in 
fourth order of J/U to C([l, 1, 0]). Because one operator 
Hi = ajaj appears in each process chain, the required 
order of perturbation theory is n = + 1, where v sig- 
nals the number of ordinary tunneling processes i?tun, as 
before. When determining the weight factor of a given 
diagram of this kind, no division by the number of sites 
occurs, because i and j distinguish specific sites. 

We have concentrated our investigations on correla- 
tions along a line parallel to a principal axis of the lat- 
tice, and along a diagonal in a main lattice plane, as 
corresponding to lattice vectors 

= [s, 0, 0] , s = l,2,...,6, 

= [t, t, 0] , i= 1,2,3, (32) 

where s and t are given as multiples of the lattice con- 
stant. Depending on j, only even or only odd orders 
contribute. Our current implementation is able to han- 
dle the expansion up to 10th order in J/U for dimen- 
sionalities d = 3 and d — 2, and up to 11th order for 
d = I. For the 3D system the number of the diagrams 
encountered is stated in Appendix |^ As in the case 
of the energy correction, the coefficients grow approxi- 
mately exponentially with the number u of ordinary tun- 
neling processes. Our findings for the ID system with 
unit filling (g = 1) again perfectly match the expansion 
reported in Ref. [20- Moreover, once again the data can 
be scaled such that they become almost independent of 
the filling factor g. To this end, we divide the corre- 
lation function by the leading density dependence, and 
plot C{fij) = C{fij)/ y/g{g + 1) for fixed scaled tunnel- 
ing parameter yJg{ci+V)J /U vs. j = jr^ j|, employ- 
ing the Euclidean norm. At least for sufficiently small 




diagonal 

1 2 3 4 5 6 

Distance r 

FIG. 8: (Color online) Logarithm of the scaled atom-atom 
correlation function C{r) = C(r)/ \/ g{g + 1) for d = 3 at 
J/U = 0.01/ g{g + 1), with various filling factors g. The 
corresponding data for d = 2 and J/U — Q.Q2/ ^ g[g + 1) 
are shown in the inset. Due to the scaling, data points for 
different g lie almost on top of each other. The decay of the 
correlations is quite well described by exponential functions, 
as testified by the linear fits. Correlations along the diagonal 
(dotted lines) decay quicker than those parallel to a main axis 
(full lines). 



^yg{g+l)J/U, we then find a beautiful exponential de- 
cay of the correlations with distance, depending on the 
direction considered: 

C(f,j) cx exp( - a{J/U)r,,,) . (33) 

In Fig. [S] we display logarithms of such scaled atom-atom 
correlations C{rij) for d = 3 and d — 2, together with 
linear fits. The correlations along the lattice axis are 
slightly larger than those along the diagonal. Figure [S] 
depicts C{fij) for the ID case, for three scaled tunnel- 
ing parameters. As expected, lower tunneling rates lead 
to quicker decays of the correlations. Again the scaling 
works remarkably well here, mapping the data for differ- 
ent filling factors nearly onto each other. Such an expo- 
nential decay of ID correlations in the regime of low tun- 
neling rates has also been observed with DMRG methods 
for distances up to 20 lattice constants by KoUath et ai— 

Figure [TO] shows the evolution of the two decay con- 
stants a (parallel and diagonal) with increasing J/U for 
the 3D and the 2D system with unit filling, determined 
from the slopes of linear fits to logarithmic plots simi- 
lar to Fig. [51 For large J/U the data should be con- 
sidered as tentative only, since the quality of the fit de- 
teriorates then. While the perturbative expansion can- 
not be expected to be valid beyond the critical hopping 
strength which marks the transition to a superfluid (with 
g = 1, one finds {J/U)c ~ 0.034 for the 3D syst em and 
'( J/f7)c « 0.059 for the 2D case, see Refs. MlM^ and 
Sec. |V|, and the true decay constants are supposed to 
vanish at that point, it is interesting to observe in Fig.fTOl 
that the tentative data obtained for the two directions 
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J/U = 0.05/(g(g+1))"2 
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"i... 



J/U = 0.15/(g(g+1)) 



1/2 



1 2 3 4 5 6 

Distance r 

FIG. 9: (Color online) Scaled atom-atom correlation C{r) — 
C{r)/^g{g+ 1) for the ID system with g = 1,2,3,5,50 at 
J/U 



0.05/^ gig + 1) (full line), J/U 
(dashed line), and J/U 



0.10/Vff(<7 + 1) 
0.15/^5(5+1) (dotted line). The 
lines are exponential fits of the form /3exp(— ar), with pa- 
rameters a and /3 determined for g = 4. With increasing 
tunneling parameter J/U the quality of the fit becomes less 
good. 



3 



o 




2D - 




+ Q 
















3 


0.025 


0.05 



parallel + 
diagonal o 



0.01 0.02 
J/U 
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FIG. 10: (Color online) Correlation decay constants a{J/U) 
for d = 3 vs. tunneling coupling J/U, as determined tenta- 
tively from fits to the correlation functions. Observe that 
the coefficients for different directions (parallel and diagonal) 
converge with increasing {J/U). The inset shows the corre- 
sponding results for d = 2. All data have been computed for 
unit filling [g — 1). 



converge with increasing J/U. 



C. Density-density correlations 

Similar to the atom-atom correlation Cifij), we inves- 
tigate the density-density correlation 

Dir.^j) = = {n^nJ) . (34) 

Besides the nearest-neighbor tunneling processes iftun, 
every process chain now contains one operator Hi — 
fiifij, sketched in the diagrams by two diamonds con- 
nected by a line, as illustrated in Fig. [TT] Because this 



^ 



a) 



b) 



c) 



FIG. 11: Diagrams for calculating the density-density cor- 
relation D(r). Subfigures a) and b) depict the second-order 
diagrams required for D([l, 0, 0]), with associated weight fac- 
tors a) 2(2d — 1) and b) 1. The linked diamonds denote the 
operator hihj, while arrows represent nearest-neighbor tun- 
neling. Subfigure c) depicts a more complicated diagram for 
D([l, 1, 0]) of 6th order in J/U. 



process Hi does not change the particle number and 
leaves the state it acts on unaltered, the diagrams for 
Di j always contain an even number of ordinary tunnel- 
ing processes, as necessary for generating closed loops, 
and therefore only even orders in J/U contribute. Since 
(m|_ffi|m) does not vanish, we cannot reduce the num- 
ber of Kato terms as much as was possible in the case 
of the energy correction and the atom-atom correlation, 
leaving us with higher computational effort. Moreover, 
the number of diagrams is larger than in the previous 
situations, as shown in Appendix \K\ for d — 2>. We com- 
pute density-density correlations up to order eight in the 
tunneling parameter J/U for the 2D and the 3D system, 
and up to order ten in the ID case. 

For d = 3 the corrections AD{fij) = D{fij) — g"^ 
to the zeroth-order value 5^ obtained for J/U ~ are 
fairly small, as exemplified in Fig. [T^l In the BEC-limit 
[J/U — > cxd), the density-density correlations again are 
given by Dij — g^, assuming large systems {M ~+ 00). 
For small tunneling parameter J/U we observe an expo- 
nential decay of AD{rij) with increasing distance r^j , as 
in Fig. [T^l Similar to the case of the atom-atom correla- 
tions, the decay constants depend on the direction: "Di- 
agonal" correlations tend to decay faster with distance 
than "parallel" ones. 



D. Occupation number fluctuations 

The squared fluctuations of the site-occupation num- 
bers are given by the variance 



{Ahf 



(35) 



Due to the homogeneity of the Bose-Hubbard system, 
this quantity is independent of the site index, and the 
expectation value of the number operator n is just the fill- 
ing factor g. Thus, for determining the variance (j35p we 
need to know (n^), and therefore generate our diagrams 
such that besides ordinary tunneling processes i?tun one 
operator Hi = nf appears. Because this process does 
not alter the system's state, again all Kato terms have to 
be evaluated, as in the case of the density-density corre- 
lation. Although the diagrams now look very similar to 
the ones encoding the energy correction, their number is 



9 





- 1 U 














(N 

cn 


-20 




+'■■••.. 






g=2 


X 






g=5 


■ 




-25 


■ g=io 


Q 






g=20 


• 






g=50 






-30 


g=100 


A 






parallel - 








diagonal ■ 





-35 



-10 
-20 





2D 











"a. 



Distance r 



FIG. 12: (Color online) Logarithm ol the correction AD{r) — 
D{r) — to the zeroth-order density-density correlation for 
d = 3 at J/U = 0.01/ \/g(g + 1). Because of this scaling, data 
points for different filling factors g lie almost on top of each 
other. The decay of these corrections is quite well described 
by exponential functions, as emphasized by the l inear fits . 
The inset shows data for d = 2 with J/U = 0.02/ ^g{g + 1). 
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FIG. 14: (Color online) Occupation number fiuctuations An 
for filling factors g = 1, 2, 3, 5, 10, 20, 50, and 100 for the 
3D (lower data points), the 2D (upper data points), and the 
ID system (inset). Vertical lines mark the critical hopping 
strength {J/U)c for the Mott insulator-to-superfiuid tran- 
sition with unit filling. Plotted vs. the scaled parameter 
2dy^g(g + 1)J/U, data points for different filling factors fall 
onto each other. For low hopping strength, the fluctuations 
grow linearly with J/U. 



THE MOTT-SUPERFLUID PHASE 
TRANSITION 



FIG. 13: Set of diagrams of sixth order in J/U for calculating 
(nf). The operator is marked by a diamond, which is 
added here in four topologically different ways to a sixth- 
order diagram for the energy correction. This leads to four 
difi'erent diagrams contributing to the perturbation series in 
seventh order. 



much higher when considering equal numbers v of tun- 
neling processes, as revealed by Tab. |V] in Appendix [X] 
The reason for this increase is evident in Fig. [T31 The 
topology of the diagrams becomes more complex by in- 
troducing the additional operator n|, depicted by a di- 
amond at site i. In the example shown in Fig. 1131 one 
diagram contributing in sixth order of the tunneling pa- 
rameter J/U to the energy correction gives rise to four 
different diagrams for the calculation of (n^). 

We were able to determine the expansion for (tt?) up 
to order ten in the tunneling parameter J/U for dimen- 
sionalities d— 1, 2, 3. The number fluctuation An grows 
approximately linearly with J/U for J/U < {J/U)c, and 
our scaling for different filling factors once more works 
very well in this parameter regime, as demonstrated in 
Fig. [TH Since the critical value {J/U)c for the Mott 
insulator-to-superfluid transition is roughly proportional 
to 1/g, this implies that the relative fluctuation An/{h) 
behaves like 1/g close to the transition point. At {J/U)c 
we find An « l/2d. In the superfluid regime, where the 
expansion in J/U is no longer valid, the fluctuation An 
eventually approaches the value ^/g. 



A further fruitful application of the diagrammatic 
many-body perturbation theory based on Kato's se- 
ries ([5]) consists in the accurate determination of the 
boundary between the Mott phase and the super- 
fluid phase for the homogeneous Bose-Hubbard modelj^ 
Qualitatively, the phase diagram of the Bose-Hubbard 
model has been understood since the late eightiesi^ More 
quantitatively, it has been intensely studied, e.g., by 
means of the strong-coupling expansion conducted by 
Freericks, Elstner, and Monien^i^ in one, two, and three 
dimensions. Recently, the quantum Monte Carlo analy- 
sis by Capogrosso-Sansone et ali^^^ has provided quasi- 
exact values for g = \. In one dimension, fairly large 
systems even including a confining trap potential can be 
treated with DMRG techniques>i^iiiii^ii^ But so far, es- 
pecially for dimensionalities d > 1 it has remained hard 
to obtain precise results for filling factors well above 
g ^ 1. Our approach is able to fill this gap. 



A. Method of effective potential 

For locating the parameters {J/U)c marking the quan- 
tum phase transition, we employ the method of the effec- 
tive potential,— in the formulation recently given by dos 
Santos and Pelsteri^ To begin with, one adds spatially 
constant source and drain terms to the Bose-Hubbard 
Hamiltonian ([2]), such that particles are created and an- 
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nihilated with uniform strengths r] and rj* at each site: 



(36) 



Since our considerations apply for any fixed value of the 
chemical potential, we do not explicitly indicate the de- 
pendence on /i/J7 in the following. We now define the 
grand canonical free energy at zero temperature 



F{J/U,rj,7j*) ^ {Hbb) 



(37) 



as the ground-state expectation value of the full Hamil- 
tonian pop for finite source strength, and expand this 
expression into a power series in rj and r/*: 



F{J/U,ij,7j*) - M fo{J/U) + y^C2n{J/UM 



2n 



71 = 1 



(38) 

The appearance of only powers of \ri\'^ reflects the fact 
that the free energy does not depend on the phases of 
and rj* . The intensive quantity /o denotes the free energy 
per lattice site in the absence of the sources. The coef- 
ficients appearing in the above expansion, in their turn, 
are expanded into power series in the hopping strength 
J/U, giving 



C2nij/u) = y2^2rhj/ur ■ 



(39) 



Whether the system is a Mott insulator or a superfluid 
is determined by its reaction to the sources. Hence, we 
introduce the functions 



(40) 



where the respective second equalities are consequences 
of the Hellmann-Feynman theorem. Assuming the in- 
vertibility of these functions, we then perform a Legen- 
dre transformation from _F to a function F depending on 
J /U , ip, and 'ip* as independent variables: 

r{J/U,ilj,ip*) = F/M ~^j*r]~iPr]* . (41) 
From Eqs. and ((55)) one obtains 

V'(?7,r?*) ^C2V +2c4|?7pr/ +0(77^) , 
riv, V*) =C2V* + 2c4|?7|\*+0(?75) . (42) 

Inverting these relations, and inserting into Eq. (HII), one 
arrives at an expansion of F in powers of \'4'\'^' 

T{J/U,i^,r) = fo - + ^1^1^ + 0(|^n . (43) 

Because rj and ip* , as well as rj* and ip^ constitute Leg- 
endre pairs, one also has the identities 



or 

'Wp* 



-rj and 



(44) 



Now the original Bose-Hubbard system ([2]) is recovered 
from the extended system ([55)) by setting rj = ij* = 0. 
Hence Eq. ()44p implies that the system adopts that value 
-00 which renders F stationary. This is akin to a mechani- 
cal system adopting a configuration in which its potential 
is stationary, signaling the absence of external forces, and 
thus motivates to dub F as an "effective potential" . 

Unless fi/U is integer, one finds C2 < for sufficiently 
small J/U, whereas C4 > (see Appendix IB)) . so that 
one has ^pQ — 0; this is characteristic for the Mott phase. 
Upon increasing J/U, the order parameter ^0 takes on 
a non-zero value when the system enters the superfluid 
phase, indicating long-range phase coherence. Hence, 
for any given value of the chemical potential the phase 
boundary (J/t/)pb is determined by that J/U for which 
the minimum of the expression (|43p starts to deviate from 
iV'oP = 0. Evidently, this occurs when the coefficient 
~l/c2 of IV'P vanishes. We point out that C2 can be re- 
garded as a susceptibility^^ Xi being the derivative of the 
function ip{ri, rj*) with respect to the source 77: 



dtp 



C2 • 



(45) 



In effect, one has to identify that hopping parameter J/U 
for which the susceptibility C2 diverges; this divergence 
marks the quantum phase transition. 

In order to compute C2 by means of the process chain 
approach, we add the perturbation 



v^-j/uy: 



a] a. 



E 



(46) 



to the zeroth-order Hamiltonian i/o, as implied by the 
extended system (1551) . Since C2 is the coefficient of jyyp 
in Eq. (|55)) . it follows by comparison of coefficients with 
Kato's series (|8]) that only chains containing one creation 
process {rja\) and one annihilation process {rj*(Xj) con- 
tribute to C2. We adjust our diagrams by introducing 
a creation process, symbolized by a dot (•), and an an- 
nihilation process, indicated by a cross (x). Since the 
operations of creation and annihilation alter the parti- 
cle number, the tunneling processes do not need to form 
closed loops here, in contrast to the cases examined be- 
fore. This leads to contributions in even and odd orders 
of J/U . For constructing the diagrams with a specified 
number v of tunneling processes, and for appending the 
correct weight factors, we generate all paths from an ini- 
tial site to any other site which can be reached with v 
nearest-neighbor tunneling events. The number of such 
paths behaves like (2^)", which is to a good approxima- 
tion equal to the sum of all weight factors. As examples 
for the emerging diagrams, the first orders = 0, 1, 2, 
and 3 in the tunneling parameter J/U are visualized in 
Fig.im 

Table [Hi which lists the number of diagrams for the 
2D and the 3D system, shows that these numbers re- 
main equal for both cases up to order v = 1 (the weight 
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FIG. 15: Diagrams for determining C2 up to order 3 in the 
tunneling parameter J/U. Subfigure a) shows the zeroth- 
order diagram with creation (•) and annihilation (x) taking 
place at the same lattice site, and weight factor 1. In b) 
we depict the only first-order diagram; its weight factor is 2d. 
The second-order diagrams in subfigure c) have weight factors 
2d and 2d{2d — 1). The weights of the third-order diagrams 
in d) are (from left to right) 2d, 2d{2d - 1), 2d{2d - 1), and 
2d{2d-lf. The one-way diagrams, which acquire the largest 
weights for high dimensionality, are diagram a) , b) , the second 
diagram in c), and the last diagram in d). 



TABLE II; Number of diagrams to be evaluated when cal- 
culating the phase boundary for the 2D and the 3D Bose- 
Hubbard model to i^th order in the hopping parameter J/U, 
corresponding to the order 1^ + 2 of Kato's perturbation series. 



V 
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10 


d = 
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4 


10 


22 


58 


140 


390 


988 


2815 


d = 


3 


1 


1 


2 


4 


10 


22 


58 


140 


394 







factors, of course, do depend on the lattice dimension- 
ality). The first difference occurs in eighth order in the 
tunneling parameter, because with eight tunneling pro- 
cesses it becomes possible to construct diagrams which 
connect more than four nearest neighbors on a hypercu- 
bic lattice. By analogy, the first difference in the number 
of diagrams between the 3D and the 4D model occurs for 
1/ = 12. 



B. Mean-field limit 



For high dimensionality d the one-way diagrams, which 
avoid "back and forth tunneling" (see Fig. [T5)) . dominate 
the contributions, because their weight factors go with 
{2dY to leading power of d, whereas all other diagrams 
contribute with lower powers. Therefore, in the limit of 
large dimensionality in each order v only the one-way di- 
agram has to be taken into account, all others possessing 
negligible weight factors then. These one-way diagrams 
can easily be evaluated analytically in every given order. 
Since they are one-particle reducible, they factorize into 



their one-particle irreducible parts as followsi ^^'^" 
.X = (-l)°(.x)i 

X = (-l)^(.x)^ (47) 



Identifying the zeroth-order term ("x) with 03*^^ and 
accounting for the factor 2d which counts the possible 
directions on a d-dimensional hypercubic lattice, the rep- 
resentation of the susceptibility takes the form 



C2 



(48) 



Because this series is geometric, its radius of convergence, 
and hence the phase boundary, can be immediately read 
off from the relation 



2da 



(0) 



(49) 



pb 



Therefore, it only remains to compute a^2^ by evaluat- 
ing the diagram • x . This gives rise to two permutations, 
which we write as (x , •) and (• , x): Either the creation 
process precedes that of annihilation, or vice versa. The 
only relevant Kato term now is (m|l/5^T^|m); the respec- 
tive energy denominators enforced by the linking opera- 
tor S*^ are AiJparticie = -Bm^ - iJparticic = fJ-/U - g for an 

extra particle, and AE^hoic — Em — -Ehoic = ^l^-fU + g— 1 
for an extra hole. Thus, the two contributions figure as 



(x,.): y^TT 



1 



■particle 



.9 + 1 
l^/U-g' 



(•>x): V9- 



1 



Airhole 

combining them yields 
.,(0) _ 



-fi/U + g-l 



^l/U+l 



{fi/U-g+l)ig-,x/U) 



(50) 



Putting everything together, the phase boundary in the 
limit of infinite dimensionality, as determinded from the 
radius of convergence of the series (|48p , is given by 



2d 



pb 



i^l/U-g + l)ig-^i/U) 

^i/u + l 



(51) 



which agrees exactly with the result of the mean-field 
calculation by Fisher et al^ This not only clarifies why 
the mean-field limit coincides with that of infinite dimen- 
sionality, but also gives a vivid illustration of our general 
approach. 
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FIG. 16: (Color online) Logarithm of the coefficients 
obtained with filling factors g 
with g — 10 for the 2D case. 



1, 50 for the 3D model, and 



C. Results for lower dimensions 

For lattice dimensionalities d = 2 and d — 3, we have 
computed the (negative) coefficients a'^^ up to order v — 
8. Again, these coefficients grow to good approximation 
exponentially with the number v of tunneling processes 
taken into account, as illustrated in Fig. [11] for filling 
factors g — I, 50 {d — 3) and 10 [d = 2). Data for other 
filling factors behave similarly. 

In the case of infinite dimensionality, the ratio 



stays constant, and directly yields the phase 
boundary. In contrast, for finite d this ratio still changes 
slightly with increasing ly. We therefore make use of a 
simple extrapolation scheme, based on d'Alembert's ra- 
tio test:'^^ The radius r of convergence of a power series 
s — J2^=o b'-'^^'z" is given by 



lim 



(52) 



if this limit exists. Therefore, we determine the phase 
boundary (J/C/)pb by plotting the ratios l^'i^ fo^' 

orders = 1 to 8 vs. l/i^, and by extrapolating to v ^ 
oo by means of a linear fit. Figure [T7] illustrates this 
scheme for dimensionalities c? = 2, 3, 5, and 10, assuming 
unit filling. Note that the slope of the straight fitting 
lines decreases with increasing dimensionality, signaling 
the approach to the strictly geometric series present for 
d = oo. 

This method of extrapolation also provides a reliable 
estimate of the systematic error. If we employ different 
selections of coefficients aj"'' , such as those with i/ = 2 to 
8 or 1/ = 3 to 8, we obtain very similar values for the phase 
boundary, thus confirming the high fidelity of our data. 
This is shown in Tab. IIIIl which lists raw data for the 
critical hopping strengths {J/U)c, marking the position 
of the tip of the respective Mott lobe. For the 2D system 
we thus estimate the overall relative error to be smaller 
than 2%, while it is smaller than 1% for the 3D model, 
and reduced still further in the 4D case. The compari- 
son of our results with recent data for g = 1 obtained by 
QMC method o '^-' — shows a remarkable agreement. Only 
close to the tip of the lobes small deviations are visible. 
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FIG. 17: (Color online) Extrapolation scheme for determining 
the phase boundary: The ratios a'^'^'^ /a'^^ are plotted vs. 
l/iy, and extrapolated linearly to 1/;/ = 0. Data are given for 
dimensionalities d = 2, 3, 5, and 10, and chemical potential 
fi/U — 0.5. Observe that the slope decreases with increasing 
dimensionality. 



TABLE III: Critical values of the hopping parameter {J/U)c 
for dimensionalities d = 2, 3, and 4, and filling factor g = 1. 
These data were obtained from linear fits to coefficients from 
different orders v, as stated in the left column. 









d=2 


d = 3 


d = 4 


1 - 8 


5.9093E-002 


3.4068E-002 


2.4131E-002 


2 - 8 


5.9853E-002 


3.4248E-002 


2.4189E-002 


3 - 8 


5.9403E-002 


3.4092E-002 


2.4107E-002 


4 - 8 


5.9846E-002 


3.4255E-002 


2.4163E-002 


5 - 8 


5.9482E-002 


3.4080E-002 


2.4093E-002 



as revealed by the inset of Fig. [18] for d = 3. A similar 
comparison for d = 2 can be found in Ref. [H. The higher 
the filling factor, the more pronounced the model's ap- 
proximate particle-hole symmetry-S. becomes, which ren- 
ders the Mott lobes symmetric with respect to the chem- 
ical potential, such that the critical chemical potential 
approaches (/x/J7)c = g — 0.5 for high g (see Fig. [TSt . 
Interestingly, when multiplying the critical hopping pa- 
rameters {J /U)c for fixed dimensionality d and varying g^ 
as listed in Tab. HVl by g{g + 1), all data fall within a 
quite narrow range, as witnessed by Fig. 1191 

We point out that the calculation of the phase bound- 
ary for d = 1 requires further considerations, because of 
a re-entrance phenomenon>i£ For certain values of the 
chemical potential the transition from the Mott insula- 
tor to the superfluid in the ID system is followed by an- 
other transition back to the insulator phase upon increas- 
ing J/U, before the superfluid phase is reached again. 
Thus, for one value of n/U there then exist three phase- 
bounding values of J/U, which cannot immediately be 
extracted with our present procedure. 



D. Approaching the mean-field limit 

Although of lesser experimental relevance, it is still 
interesting to investigate systems with dimensionality 
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FIG. 18: (Color online) Mott lobes for d = 3, and various 
filling factors g. The dashed line marks the limit [fi/U)c = 
g — 0.5 attained for high g. For g = 1, QMC data^^ are 
included. The inset magnifies the tip of the lobe belonging to 
unit filling. 



TABLE IV: Critical values {^^/U)c and {J/U)c for various 
filling factors g. For locating the tip of the respective Mott 
lobe, /x/C/ has been varied in steps of 0.001. Relative errors 
of {J/U)c are less than 1% for d = 3, and less than 2% for 
d = 2. 





d = 2 


d = 3 


9 






(M/f/)c 


{J/U)c 


1 


0.376 


5.909E-002 


0.393 


3.407E-002 


2 


1.427 


3.480E-002 


1.437 


2.007E-002 


3 


2.448 


2.473E-002 


2.455 


1.427E-002 


4 


3.460 


1.920E-002 


3.465 


1.108E-002 


5 


4.470 


1.569E-002 


4.472 


9.055E-003 


6 


5.472 


1.327E-002 


5.476 


7.657E-003 


7 


6.476 


1.150E-002 


6.479 


6.634E-003 


8 


7.479 


1.014E-002 


7.482 


5.852E-003 


9 


8.481 


9.073E-003 


8.484 


5.235E-003 


10 


9.483 


8.208E-003 


9.485 


4.736E-003 


20 


19.491 


4.202E-003 


19.492 


2.425E-003 


50 


49.496 


1.706E-003 


49.497 


9.842E-004 


100 


99.498 


8.571E-004 


99.498 


4.946E-004 


1000 


999.500 


8.609E-005 


999.500 


4.968E-005 


10000 


9999.500 


8.613E-006 


9999.500 


4.970E-006 



d > 3, in order to study the convergence towards the 
mean- field limit. For higher d it becomes harder to ob- 
tain the necessary diagrams, and their weight factors. 
Nonetheless, we are able to treat systems of arbitrary 
dimensionality at least up to order v = A in the tun- 
neling parameter J /U , because the corresponding weight 
factors can still be figured out by combinatorial reason- 
ing. Despite this relatively low order the precision of the 
phase boundaries thus obtained is quite high for large d, 



since the fluctuation of the ratio ^'/q;2'''' decreases 
significantly with increasing d, as illustrated by Fig. 1171 
The resulting Mott lobes for unit filling are displayed in 
Fig. [20] for d = 2, 3, 5, 10, and 20, together with the 
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FIG. 19: (Color online) Scaled values ^ g{g + l){J/U)c of the 
critical hopping parameter vs. filling factor g for d = 2 (upper 
panel) and d = 3 (lower panel). Observe the rather fine scale 
at the left margin. 
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FIG. 20: (Color online) Phase diagram for dimensionalities 
d = 2, 3, 5, 10, and 20, together with the mean-field phase 
boundary, for g — 1. With increasing dimensionality the 
data approach the mean- field prediction (|5ip . In the inset the 



scaled critical hopping strength 2dy^g{g + l){J/U)c is plotted 
vs. 1/d for filling factors g = 1 (lower data set) and g — 10000 
(upper data set). 



mean-field phase boundary. In the limit d ^ oo the crit- 
ical parameter {J/U)'^^ can be deduced from Eq. (|5ip . 
giving 



2d 



mf 



2g+l-2,/g{g+l) 



(53) 



which scales like 1/g for large g. Hence, the approach 
to the mean-field limit can be well monitored by plot- 
ting the phase-bounding chemical potentials for each d 
vs. 2dJ/U, as in the main part of Fig. [20l While the 
curves agree fairly well with each other at the edges of 
the lobes even for low dimensionalities, larger deviations 
occur around the tips. The convergence to the mean- 
field phase boundary with increasing d is clearly visible. 
When plotting the scaled data 2dyJ g{g + l){J/U)c for 
fixed g as functions of 1/c? in order to directly highlight 
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the approach to infinite dimensionahty, as done in the in- 
set of Fig. [^nifor 5 = 1 and g — 10000, we obtain smooth 
curves aiming for 1/c? ^ at the respective mean-field 
result determined by Eq. ([55]) . The scaled data belonging 
to different g fall into a remarkably narrow range; their 
relative separation amounts to about 3%. In fact, for 
any dimensionality d and any filling factor g the critical 
values computed in this work are fitted by 



mf 



0.13 



V.9(5 + 1) ' 



(54) 



with an accuracy of about 1%, where (J/C/)™^ follows 
from Eq. ([53]) . Note that even for relatively high d the 
deviation from the mean-field prediction is not negligible. 
For example, even for c? = 10 and g — 1 the value of 
{J/U)c still exceeds the mean-field limit by about 4%. 



VI. CONCLUSION & OUTLOOK 

The essence of the high-order process chain approach 
is captured by the example given in Subsec. flll Al On the 
one hand, one has to generate the Kato terms, as dictated 
by the perturbation series ([8]). This step is universal, and 
thus has to be done only once for all kinds of perturbative 
calculations. On the other hand, one has to construct all 
the diagrams pertaining to the particular observable un- 
der study, and their weight factors. These diagrams then 
are worked out in a Cinderella-type fashion: Each per- 
mutation of the processes constituting a diagram has to 
be compared to the pattern of intermediate states pro- 
vided by the Kato terms; the matching permutations are 
evaluated, the others discarded. The bottlenecks of this 
scheme are the generation of the diagrams, which poses 
nontrivial combinatorial tasks in higher orders, and the 
factorial growth of the number of permutations with the 
order of perturbation theory. While it might be feasible 
to optimize diagram-generation with the help of specifi- 
cally adapted routines, the explosive growth of the num- 
ber of permutations currently appears to limit straight- 
forward numerical applications of this algorithm to about 
the twelfth order. 

But still, this could open up substantial possibilities 
for the further analysis of strongly correlated quantum 
many-body systems. This suggestion is underlined not 
only by our calculations of the various ground-state cor- 
relation functions for the homogeneous Bose-Hubbard 
model in Sec. lIVl but also by the comprehensive determi- 
nation of its phase diagram in Sec. |Vl The entire set of 
all dimensionalities d> 2 and all filling factors g has been 
covered by a single approach, giving excellent agreement 
with previous findings in those cases for which accurate 
calculations had been performed before. 

The conceptual ease with which these results have been 
obtained here suggests that the process chain approach^^ 
should also turn out useful for the theoretical investi- 
gation of other systems which so far are less well un- 
derstood. We expect our strategy to work with similar 



success for different types of lattices, such as triangu- 
lar or hexagonal ones, for ladder systems;^ and for lat- 
tices with a superstructure, such as recently considered 
in Ref . M- 



Acknowledgments 

We thank M. Langemeyer and V. Steenhoff for im- 
portant help with the generation of high-order Kato 
terms, and B. Capogrosso-Sansone for providing the 
QMC dataii^iii Computer power was obtained from the 
GOLEM I cluster of the Universitat Oldenburg. N. T. 
acknowledges a fellowship from the Studienstiftung des 
deutschen Volkes. A. E. thanks M. Lewenstein for kind 
hospitality at ICFO-Institut de Ciencies Fotoniques and 
acknowledges a Feodor Lynen research grant from the 
Alexander von Humboldt foundation, as well as sup- 
port by the Spanish MEC (grant FIS2008-00784 "TO- 
QATA", ESF-EUROQUAM program FIS2007-29996-E 
"FERMIX"). This work was further supported by the 
Deutsche Forschungsgemeinschaft. 



APPENDIX A: NUMBER OF DIAGRAMS 



TABLE V: Number of diagrams encountered to i/th order in 
the tunneling coupling J/U when calculating the quantities 
considered in Sec. IIVI for the 3D Bose-Hubbard model. An 
obvious shorthand notation is used here for the lattice vec- 
tors introduced in Eq. (|32p . such that C(s) = C([s, 0, 0]) and 
D{t,t) = D{[t,t,0]). 



V 
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6 
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10 


E 






1 




3 




7 




29 




127 


C(l) 




1 




3 




10 




50 




281 




C(2) 






1 




3 




15 




102 




795 


C(3) 








1 




3 




18 




102 




C(4) 










1 




3 




19 




151 


C(l,l) 






1 




3 




17 




126 




1051 


C(2,2) 










1 




3 




21 




190 


D(X) 


1 




2 




8 




40 




250 






m) 


1 




1 




6 




28 




194 








1 




1 




5 




23 




144 






D(l,l) 


1 




1 




7 




32 




227 






D(2,2) 


1 




1 




5 




22 




140 






1 




1 




4 




18 




106 




697 



APPENDIX B: POSITIVITY OF a 



The method of the effective potential outlined in Sub- 
sec. IV Al cruciallv requires that the coefficient of in 



Eq. (|43)l . and hence C4, be positive. For evaluating C4 
within the process chain approach, we have to construct 
diagrams containing exactly two annihilation and two 
creation processes. To zeroth order in J/U (fourth order 
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of the perturbation series) , one gets one diagram with two 
creations and two annihilations at the same site (•• x x), 
leading to 4!/(2!)^ = 6 permutations. Similar to Eq. (fT5|) . 
the corresponding Kato terms are {m\V S'^V S^V S^V\m) 
and {in\V S^V S^V S^V\in) . The relevant energy denom- 
inators are as follows: 

One hole: AEh ==.g - 1 - fi/U 

Two holes: AEhh=2g - 3 - 2/i/C/ 

One particle: AEp ^fi/U — g 

Two particles: AEpp^2fi/U - {2g + 1) . 

All these energy differences are negative, with 



AEh > AEhh , 
AEh > AEpp , 



AE„ > AE, 



AEr, 



pp 1 
> AEhh , 



as follows from g — \ < fi/U < g. The combination of all 
contributions then yields 



5 + 1 


-g + 2 


5 + 1 


5 


{AEpf 


AEpp 


AEp 


AEh 


9 


".9-1 


5 + 1 


5 


[AEhY 


AEhh 


AEp 





(B2) 



Since according to the above relations (|B1[) both factors 
in square brackets are positive, the coefficient C4 is man- 
ifestly positive to zeroth order in the tunneling parame- 
ter J/U. We have investigated higher orders in J/U nu- 
(Bl) merically, and obtained only positive contributions a^J^\ 
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